Exact few-body results for strongly correlated quantum gases in two dimensions 
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The study of strongly correlated quantum gases in two dimensions has important ramifications 
for understanding many intriguing pheomena in solid materials, such as high- Tc superconductivity 
and the fractional quantum Hall effect. However, theoretical methods are plagued by the existence 
of significant quantum fluctuations. Here, we present two- and three-body exact solutions for both 
fermions and bosons trapped in a two-dimensional harmonic potential, with an arbitrary s-wave 
scattering length. These few-particle solutions link in a natural way to the high-temperature prop- 
erties of many-particle systems via a quantum virial expansion. As a concrete example, using the 
energy spectrum of few fermions, we calculate the second and third virial coefficients of a strongly 
interacting Fermi gas in two dimensions, and consequently investigate its high-temperature thermo- 
dynamics. Our thermodynamic results may be useful for ongoing experiments on two-dimensional 
Fermi gases. These exact results also provide an unbiased benchmark for quantum Monte Carlo 
simulations of two-dimensional Fermi gases at high temperatures. 



I. INTRODUCTION 

Two-dimensional (2D) strongly correlated quantum 
gases present unique features from the point of view of 
many-body physics [1|. Many sophisticated collective 
phenomena arise because of reduced dimensionality, such 
as the long-sought Berezinsky-Kosterlitz-Thouless transi- 
tion [2-4^1 and high- Tc superconductivity |^]. In addition, 
particles in 2D can have non-Abelian quantum statistics, 
which is strikingly different from bosons and fermions. 
For this reason, a 2D quantum system is a potential plat- 
form for topological quantum computation in a way that 
is naturally immune to decoherence [Q]. 

Recent experiments with ultracold atoms offer a unique 
opportunity to investigate this physics in a controllable 
way [1, 7\. In these experiments, one can modify aspects 
of the underlying geometry and interactions between the 
atoms, at temperatures down to one billionth of a de- 
gree above absolute zero. Experimental schemes to pro- 
duce a 2D atomic quantum gas include a one-dimensional 
(ID) optical lattice, formed by the superposition of two 
running laser waves |8l-[l2|. and strongly focused ellip- 
soidal optical traps. Using the technique of Feshbach 
resonances [13], the interatomic interaction can also be 
changed from infinitely weak to infinitely strong. This 
has already led to the observation of the crossover from 
a Bose-Einstein condensate (EEC) to a Bardeen-Cooper- 
Schrieffer (BCS) superfiuid in three dimensions [1, 7]. 

Theoretical investigations of 2D strongly correlated 
atomic quantum gases, in particular the study of super- 
fluidity in atomic Fermi gases, have already attracted 
intense attention in the past few years [ll. [l4l - |20| . How- 
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ever, theoretical methods for non-integrable 2D Fermi 
systems are limited due to significant quantum fluctua- 
tions. Although a mean-field approach combined with 
perturbation theory are usually adopted in the under- 
standing of the BCS-BEC crossover in three dimensions 
|ll [tI. l2l| - [23fl . they may simply break down in 2D. Other 
traditional methods in condensed-matter physics, such as 
exact diagonalization and quantum Monte Carlo simula- 
tion, are often less helpful than one may expect, due to 
the restriction to finite number of atoms or due to Fermi 
sign problems. Furthermore, the harmonic trapping po- 
tential in ultracold atom experiments, which is used to 
prevent the atoms from escaping, complicates theoretical 
treatments. 

In this paper, we present a few-particle perspective 
on strongly correlated 2D systems by exactly solving for 
the eigenstates of three identical fermions or bosons in 
a 2D isotropic harmonic trap, with arbitrary interac- 
tion strength. Threc-fcrmion or three-boson problems in 
three dimensions (3D) have been thoroughly investigated 
|24l - [27| , covering many aspects such as the three-body re- 
combination rate (or stability) |28l [29| , three-body per- 
spective on BEC [30], and Efimov physics [3l|, [S^l- The 
three-particle problem in low dimensions, however, is less 
well-studied despite its considerable importance. There 
are very few studies of universal low-energy properties of 
three identical bosons confined in 2D [33-36] . 

Here, by constructing the exact wave functions, we 
solve and discuss the full exact energy spectrum of three 
identical trapped fermions or bosons in 2D. As the Efi- 
mov effect occurs only when the dimensionality is greater 
than two [2J|, all the states of fermions and bosons that 
we study have universal properties determined by a single 
parameter: the s-wave scattering length Osc- For three 
bosons, we find that an attractive interaction leads to 
two distinct three-boson bound states in the form of a 
self-bound boson droplet, as predicted by Hammer and 
Son J35| using a 2D effective field theory. 



Using few-particle exact solutions, we can also solve 
the problem of a strongly correlated 2D quantum gas 
at high temperatures, includin g b oth thermodynamics 
|37| and dynamical properties [38, |39|, using a quan- 
tum virial expansion method [40]. Here, we calculate 
the second and third virial (expansion) coefficients of a 
2D Fermi gas. We then investigate the high-temperature 
equation of state, including the chemical potential, en- 
ergy and entropy, as a function of temperature at a given 
interaction strength. Our thermodynamics results give 
valuable insights for ongoing experiments on 2D Fermi 
gases[4lL l42j. Further, these results may also provide a 
useful benchmark for quantum Monte Carlo simulations 
for a 2D Fermi gas at high temperatures, where conver- 
gence checks are otherwise difficult to obtain. 

The paper is organized as follows. In the next sec- 
tion, we present exact solutions for the energy eigen- 
states of three-fermion and three-boson systems with ar- 
bitrary s-wave interaction in an isotropic 2D harmonic 
trap and discuss the resulting energy spectrum. In Sec. 
Ill, we calculate the second and third virial coefficients 
of a 2D Fermi gas, at a given temperature and inter- 
action strength. Then, in Sec. IV, we investigate the 
high-temperature thermodynamics of a strongly corre- 
lated 2D Fermi gas. Sec. V is devoted to conclusions 
and final remarks. In the Appendix, we outline some 
numerical details of the exact solutions. 



II. EXACT FEW-PARTICLE SOLUTIONS IN A 
2D HARMONIC TRAP 



tion strength (i.e., unitarity limit), the scattering length 
flsc in 2D is always positive due to the existence of a 2D 
bound state. 

Following the idea introduced into two-body physics 
by Bethe and Peierls [4J|, it is convenient to replace the 
s-wave interaction by a set of boundary conditions, which 
in 2D take the form [H, [M S Ell , 
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when particles i and j are close to each other. Here, 
'0(pi,--- , pn) is the wave function of a system of N 
particles and pij — \pi — pj\. In addition to these 
Bethe-Peierls boundary conditions, the wave function 
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with no two particles at the same coordinate. 



A. Two particles in a 2D harmonic trap 

As a preliminary study, let us first revisit the two- 
particle problem [43|. In a harmonic trap, the motion 
of the center-of-mass C = {pi + P2) /2 can be separated 
from the relative motion, and the relative Hamiltonian is 
given by. 



We consider a 2D few-particle system of either fermions 
or bosons in an isotropic 2D harmonic trap V{p) = 
muj^p^/2 with p = \J x^ -|- y^, where pj = {xi,yj) is the 
the j-th particle coordinate. For low-energy scattering, 
the attractive interactions between atoms can be formally 
described by a positive s-wave scattering length asc- For 
identical fermions, there is no s-wave partial wave in- 
teraction due to the Pauli exclusion principle. We thus 
consider for fermions two different hyperfine (i.e., pseudo- 
spin) states, with the interaction occurring only for two 
fermions with unlike spins. In the case of a Feshbach res- 
onance, which allows an adjustable interaction strength, 
we focus on the case of a broad rather than narrow res- 
onance. This allows us to analyse the problem without 
considering an explicit molecular channel. More gener- 
ally, the molecular field causing the resonance should be 
included, leading to a modified two-particle bound state 
eigenfunction[43] . 

A peculiar feature of 2D interactions is that any 
attraction, whatever how small, will support a two- 
particle bound state with binding energy Eb = 
4SV[exp(22)ma^J, where 7 ~ 0.577216 is the Euler 
constant [lg|. The interactions can then be alterna- 
tively characterized by the two-particle binding energy 
Eb- Contrary to the 3D BEC-BCS crossover situation, 
where the bound state appears only at a certain interac- 
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where p — pi — p2 is the relative coordinate and 
p = m/2 is the reduced mass. The energy level and 
the corresponding wave function of two-particle system 
can be written as E = Ecm + Erei and ^2p(C,/9) = 
4>2^ (C) ip2p (p) , respectively. Here, the subscript "2p" 
denotes the two-particle problem. 

The wave function of center-of-mass motion, (f)2™ (C), 
is simply the well-known wave function of 2D har- 
monic oscillators with Ecm = (^ricm + |f7Zcm| + l)huj, 
where the good quantum number ricm and rricm label, 
respectively, the number of nodes in the radial wave 
function and the angular momentum of the center-of- 
mass motion. The relative wave function should be 
solved by Hreitp2p (p) — Erei'<p2p ip) i ^^ conjuuctiou with 
the Bethe- Peierls boundary condition, \mip^Q[d/dp — 
\/{p\x\{p/asc))]'4'2p iP) ~ 0- The relative Hamiltonian 
has rotational symmetry and thus has a good quantum 
number of angular momentum mrei . Due to s-wave cou- 
pling, it is easy to see that only the nirei — branch of 
the relative wave functions is affected by the interactions, 
so we focus on this case. 

We start by considering the solutions to the free Hamil- 
tonian, without including boundary conditions. The free 
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Figure 1: (Color online). Relative energy spectrum of a two- 
particle system with rrirei = as a function of the dimen- 
sionless interaction parameter d/asc- The system goes to the 
strongly interacting limit when d/asc increases to an infinitely 
large value. 



relative Hamiltoniaii admits two types of solutions, ei- 
ther in terms of the confluent hypergeometric function 
of the first kind, exp(— p^/2d^)ii^i (— i^, l,p^/ci^), or in 
terms of the Kummer confluent hypergeometric function 
of the second kind, expi-p"^ /2tf)r{-p)U{-v,l, p"^ /(P) , 
where d = y^h/{iiuj) is the length scale of the trap, v is 
determined by Ej,ei ~ ^v + l)fiii;, and F is the gamma 
function. The first kind of Kummer function \F\ is regu- 
lar in the entire space and gives the standard wave func- 
tion of a 2D harmonic oscillator. In contrast, the second 
Kummer function U is singular at the origin. 

Now, let us include the Bethe-Peierls boundary condi- 
tion. It is easy to see that one must choose the second 
type of Kummer solution as the relative function, i.e.. 
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The parameter v or the relative energy Erei = (2i^ -I- 
Vjhijj is then uniquely determined by the boundary 
condition. Considering the property dxU {—v^l^x) — 
vU (1 — v^2,x) and the asymptotic behavior of the con- 
fluent hypergeometric function at a; — )■ 0, 
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we immediately obtain the energy equation. 
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Here, 7 ~ 0.577216 is the Euler constant and ipix) is the 
digamma function. 

In Fig. 1, we report the relative energy levels of a 
two-particle system with nirei = as a function of the 



dimensionless interaction parameter d/agc- All the en- 
ergy levels decrease with increasing interaction strength, 
as expected for an attractively interacting system. The 
lowest level corresponds to the ground state of a molecule 
with size asc and thus towards the stroirgly interacting 
limit (i.e., asc — > 0), it diverges as —h^/ima^^. All 
the other excited levels instead converge to the non- 
interacting limit. 

It is interesting to note that with a positive scatter- 
ing length, the two particles interact repulsively if they 
do not occupy the ground state of molecules. Thus, by 
excluding the lowest energy level. Fig. 1 can be alterna- 
tively viewed as the energy spectrum of two repulsively 
interacting particles [4^. Then, the right side with van- 
ishing asc is the non-interacting limit for the repulsively 
interacting system, and the unitarity limit of infinitely 
large asc is the strongly interacting limit. 

In the limiting case of either zero or infinite scattering 
length, one may calculate the asymptotic behavior of en- 
ergy levels. We find that, for the n-th level, the relative 
energy is given by, 
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where n = 0,1,2,... is a non- negative integer and, in 
the limit of asc -^ 0, the lowest molecule state has been 
excluded in the count of energy levels, so that 71 = 
correspoirds to the first excited state. 



B. Three fermions in 2D harmonic trap 

Let us irow turn to the three-particle problem. For 
three fermions, we consider the configuration with two 
spin-up fermions (particle 1 and 3) and one spin-down 
fermion (particle 2), i.e., a fit configuration. It is 
convenient to use Jacobi coordinates. We define the 
center-of-mass coordinate pcM = (pi + P2 + Ps) /3, to- 
gether with two relative coordinates r = p\—p2 and 
p = (2/V3) [p3 - (pi +p2)/2]. The solution for the 
center-of-mass motion is again the standard wave func- 
tion of a 2D harmonic oscillator. For the relative motion, 
on top of the Bethe-Peierls boundary conditions, the rel- 
ative Hamiltonian reads (26l. 
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To solve the three-fermion problem, we exteird the ap- 
proach of Efimov|3l| to the trapped case and consider 
the following relative wave function [37|. 



rsf = il-V,3)x{r,p) 



where 
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(11) 



Rnm (p) is the standard radial wave function of 2D har- 
monic oscillators with energy (2n + \m\ + l)huj, and the 
set of parameters Vm.n is determined by, 



Erei = [{2n + \m\ + 1) + (2z/™,„ + 1)] hu, 



(12) 



for a given relative energy Erei and the two good quan- 
tum numbers n and m. 

The wave function ([TU)) is easy to understand. It 
is simply a summation of products of the wave func- 
tion of the paired fermions (1 and 3), V'2p'(''; ^^m.rOj ^^^ 
of the wave function of particle 3 relative to the pair, 
Rnm (p) e^"^'^ / y/2TT . The product certainly satisfies the 
relative Hamiltonian © and gives rise to the energy con- 
servation equation for Vm.m Eq. (|12p . Owing to the rota- 
tional symmetry of the relative Hamiltonian, the angular 
momentum is well-defined and conserved. In the rela- 
tive wave function, we also include an exchange operator 
for particle 1 and 3, which ensures the symmetry of the 
wave function and ensures that the wave function satis- 
fies the Pauli exclusion principle. The set of coefficients 
a(^ can be uniquely determined using the Bethe-Peierls 
boundary conditions. We note that because of the ex- 
change operator, the two boundary conditions reduce to 
just one, since the other is satisfied automatically. 

We now examine the Bethe-Peierls boundary condition 
which should lead to a secular equation for the energy 
levels (Erei) and wave functions (a^^). Let us consider the 
first term, \inir^or{d/dr){l — 'Pi3)x{^, p)- Recall that 
Pi3X (r, p) = X (r/2 - \/3/?/2, -\/3r/2 - p/2), which is 
regular at origin. Therefore, we find. 
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On the other hand, in the limit of r ^- 0, 

r^f _ X (r, P) - X (~%/3p72, ^p/2) 
In {r/asc) In (r/asc) 
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(15) 



where effectively x {r,P)r^o = E„(-2)[7 + V'(-'^m,«) + 
ln(r/rf)]a,{i?„„(p)e*™'^/V2^. By substituting Eqs. dH]) 
and ([T5|) into the Bethe-Peierls boundary condition, it is 
easy to show that. 
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The above equation can be solved by projecting the 
left-hand side of the equation onto the expansion basis 



Rn'm {p), which is orthogonal and complete. This leads 
to the secular equation. 
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where the matrix elements are 



AI^, EE [7 -f V (-t^m.rO] ^nn' + ^7^C„„^ (19) 



and 
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(20) 
It is clear that Cnn' arises from the exchange operator 
Viz- In the absence of C„„' , the secular equation is iden- 
tical in form to Eq. ([7]), except for an additional degree 
of freedom which corresponds to the motion of particle 3 
relative to the paired fermions (particle 1 and 2). It then 
describes an un- correlated three-fermion system of a pair 
and a single particle. 

To solve the secular equation, one must imposes a 
cut-off rimax for the number of expansion functions of 
Rnm (p) ■ The accuracy of the numerical calculations can 
be improved by increasing rimax- The relative energy level 
Erei is then implicit in the secular equation via Vm.n- In 
practice, for a given relative energy level Erei, we diago- 
nalize the matrix A/ — {^„„'} to obtain all the possible 
interaction strengths d/age that correspond to this rel- 
ative energy. We then invert the relations asc{Erei) to 
calculate the desired energy spectrum (levels) as a func- 
tion of the interacting strength d/usc- The main numer- 
ical effort is to calculate the matrix elements Cnn'- We 
outline the details of this procedure in the Appendix. 

We note that, in both the two and three body cases, 
there are non-interacting solutions to the point-contact 
interaction Hamiltonian. There are many functions that 
vanish when two particles are at the same point. For 
the two-particle case, these are the ?ti > states. For the 
three-particle case, the situation is more complicated. An 
example as pointed out by Werner and Castin[25], is the 
Laughlin state: 



V'^e-^-i'--/'^' n [{^n+iyn)-{Xm+iym)]^''^ 

l<n<m<3 

These states are not included in our interacting solutions. 

C. Three bosons in 2D harmonic trap 

For three bosons we can construct a similar relative 
wave function to Eq. (|10|) . This takes the form. 



V'3t = (l + ^13+^23)x(r,p), 



(21) 



where 
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Compared with the fermion case, the only difference 
in the relative wave function is that we need to in- 
clude two exchange operators with positive sign to en- 
force the proper symmetry of the bosonic wave function 
|25| . This modifies the Bethe-Peierls boundary condi- 
tion and hence the secular equation. Otherwise, we fol- 
low the same derivation as in the fermion case. By us- 
ing r23X ir,p)=X (r/2 + v^p/2, V3r/2 - p/2) , we find 
that the secular matrix Aj, = {A^^,} takes the form. 



A^„, = [7 + ^ {->^m,n)] Snn' + (-1)™+' C„ 



(23) 



which has the same structure as v4^„/. The difference 
is that due to the additional exchange operator and dif- 
ferent sign before operators. The prefactor in the C„„' 
terms is (—1)™^ , instead of (— 1)™ /2 as in Eq. (fTOj). 

It is of importance that in two dimensions the three- 
particle bosonic wave functions we have constructed are 
universal, in the sense that all the three-boson proper- 
ties are determined by the single two-body scattering 
length [36]. This is contrary to the case in three dimen- 
sions where even in the zero-range-interaction limit, the 
Thomas and Efimov effect [3l| , results in a set of univer- 
sal three-boson bound states which are described by an 
additional three-body regularization parameter [31]. 

The absence of an Efimov phenomenon, however, does 
not imply the absence of three-body bound states. In 
free space, exactly two three-boson bound states ap- 
pear in two dimensions with an arbitrary two-body s- 
wave scattering length, in the form of boson droplets 
|3a |. The ground bound state has a binding energy 

E^g^ = 16.522688(l)£'s, while one excited bound state 

has E'-gl = 1.2704091(1)£;b. Here, Eb is the two-particle 



binding energy discussed earlier. 



D. Energy spectrum 

We now discuss the resulting energy spectrum of three 
fermions or three bosons. Typically, we set a cut-off 
^max = 128 for the number of radial wave functions 
Rnm (p) kept in the calculation. By doubling and halv- 
ing the value of nmax, we have checked that the relative 
accuracy of energy levels is less than < 10~^, except for 
the m = subspace for bosons, where the appearance of 
two three-boson bound states significantly decreases the 
numerical accuracy. 



1. Three-fermion spectrum 

Fig. 2 gives the relative energy spectrum of a three- 
fermion system at different relative angular momentum 
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Figure 2: (Color online). Relative energy spectrum of three 
trapped interacting fermions in 2D, as a function of the di- 
mensionless interaction parameter \n{d/asc)- We show the 
spectrum in different subspaces of relative angular momen- 
tum m. The ground state energy level in the subspace m = 1 
has been highlighted by a thick line. 



?7i , as a function of the interaction strength, \n{d/asc)- 
The ground state is in the subspace to = 1 due to 
the Pauli exclusion principle which prohibits all three 
fermions from interacting when tti = 0, as highlighted 
by a thick solid line. Compared with the two-body rel- 
ative energy spectrum, the energy levels are much more 
complicated. We observe two distinct energy levels with 
decreasing scattering length and therefore increasingly 
attractive interaction strengths. Some diverge to — oo 
as ajj?, while the others saturate to the limiting values 
that correspond to the non-interacting energy spectrum. 
This essential feature exactly resembles what we observed 
for the two-body relative energy spectrum shown in Fig. 
1, where the ground state of two particles diverges to 
infinitely negative energy, while the other excited states 
converge to the ideal, non-interacting spectrum. We note 
that the same feature has also been observed very re- 
cently in calculations of a trapped three-fermion system 
in 3D m. 



We may therefore identify the diverging energy level 
as the state that contains a tightly bounded pair or 
molecule, together with a fermion rotating around the 
molecule. The energy spacing of this "molecule and 
atom" state is roughly 2fvjj, accounting for the rota- 
tional degree of freedom of the fermion. Accordingly, the 
other saturating energy level is a state of three individ- 
ual fermions, which therefore should interact repulsively. 
In analogy to the two-particle case, we may regard these 
"individual atom" states as the energy states of three re- 
pulsively interacting fermions with the same (positive) 
s-wave scattering length, although there are necessar- 
ily many avoided-crossings between the "molecule and 
atom" states and the "individual atom" states. These 
appear particularly when the scattering length age be- 
comes comparable with the characteristic length scale of 
the harmonic trap, d. 

With this classification of energy levels in mind, the 
spectrum at the limiting cases of asc — > oo and asc — > 
are easy to interpret. The former is simply the en- 
ergy spectrum of three weakly attractively interacting 
fermions, which, analogous to the two-particle case, de- 
crease linearly as 1/ In(asc) with decreasing asc- The lat- 
ter, excluding the "molecule and atom" states, is the spec- 
trum of three weakly repulsively interacting fermions, 
increasing linearly as l/ln(a,5c) with increasing asc- It 
is readily seen that in these two limiting cases the en- 
ergy levels, together with their degeneracy, are con- 
nected smoothly with the spectrum of three ideal, non- 
interacting fermions. 
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2. Three-boson spectrum 

Fig. 3 presents the evolution of the relative energy 
spectrum of three bosons with increasingly attractive in- 
teraction strength. In this case, without the restriction 
of the Pauli exclusion principle, the ground state is in the 
subspace of zero relative angular momentum, to = 0. We 
highlight this again by using a thick line. The essential 
features of the spectrum are the same as in the spec- 
trum for three fermions. We observe both the "molecule 
and atom" branch and the horizontal "individual atom" 
branch, together with some avoided crossings between 
them. The latter branch may be viewed as the spectrum 
of three repulsively interacting bosons. 

However, there is an important difference, occurring 
in ground state subspace with tti = 0. The lowest two 
states in the "molecule and atom" branch are three-boson 
bound states. One is the ground state and the other 
is the lowest excited state. Their energy is lower than 
the total energy of two attractively interacting bosons 
and a third free-moving boson. In particular, the ground 
state energy is significantly lower in magnitude than the 
two-body binding energy Eb - As a result of these three- 
particle bound states, high numerical accuracy is difficult 
to obtain. As shown in Fig. 3a, the energy levels of 
the two bound states do not converge well even for the 



Figure 3: (Color online). Relative energy spectrum of three 
trapped interacting bosons in 2D at difference subspace, as a 
function of the dimensionless interaction parameter ln(d/asc). 
The ground state energy level in the subspace m = has 
been highlighted by a thick line. The numerical accuracy 
with m = is greatly suppressed due to the existence of 
the self-bound droplet-like states. We thus plot the spectrum 
with Wmax = 64 (solid lines), 128 (dashed lines), and 256 (dot- 
dashed hues), to show the slow convergence with respect to 
the number of expansion basis elements rimax- 



largest expansion basis (n^ax — 256) considered in these 
calculations. 

The two bound states describe a self-bound bosonic 
droplet formed via the attractive, short-ranged two-body 
potential, resembling the well-known bright soliton of at- 
tractive bosons in ID. Contrary to the Efimov state, these 
bound states are universal and their properties are deter- 
mined entirely by the single s-wave scattering length. 

We have estimated the binding energy of the two 
bound states at ln(d/asc) = 1 by extrapolating the en- 
ergy level obtained at a finite expansion basis to nmax — 
cx). The interaction strength is chosen to minimize the 
influence of the harmonic trap so that the size of the 
bound state (^ asc) is much smaller the trapping scale 
(~ d), while at the same time to maintain the nu- 
merical result as accurate as possible. Empirically, we 



find that the binding energy scales like, i?B3(nniax) — 
Eb3 (oo) cx rimax ■ This leads to £^^3 ~ 15.1Eb and 
Eg^ ~ 1.25Eb, which are reasonably in agreement with 
the accurate binding energies in homogeneous space, 
e'^°1 = 16.522688{1)Eb and e'^I = 1.2704091(1)£;b, as 
predicted by a 2D bosonic effective field theory [35j • The 
discrepancy, particularly for the ground state binding en- 
ergy, mainly comes from our insufficient numerical accu- 
racy. 



III. VIRIAL COEFFICIENTS OF STRONGLY 
CORRELATED FERMIONS IN 2D 

The knowledge of few-particle exact solutions provides 
a useful input for investigating the high-temperature be- 
havior of a strongly correlated quantum gas, by applying 
a quantum virial expansion to the thermodynamic prop- 
erties (37| or dynamical properties |38l . |39| . Here, we 
are interested in the high-temperature equation of state 
of strongly correlated fermions, which are now being ac- 
cessed experimentally in several laboratories. 

The essential idea of the quantum virial expansion is 
that at high temperatures where the chemical potential 
/i is strongly negative, the fugacity z = exjp(fi/kBT) = 
exp(/3/i) ^ 1 is a well-defined small parameter. We can 
therefore expand the thermodynamic potential 57 of a 
quantum system in powers of the fugacity, however strong 
the interaction strength is. Quite generally, we may write 



n = -kBTQi [z + 62^2 + • • ■ + fe„z" + . 



(24) 



where 6„ is the n-ih. (virial) expansion coefficient and 
takes the following form. 



62 = (Q2-Q?/2)/Qi, 

h:i = (Q3-QiQ2 + Q?/3)/Qi 



etc. 



(25) 
(26) 



Here, Qn — rr„[exp(— H/fcsT)] is the partition function 
of a cluster that contain n particles and the trace rr„ is 
taken over all the n-particle states of a proper symmetry. 
It is clear that Qn and hence 6„ can be calculated once the 
energy spectrum of up to n-body clusters is known. All 
the other thermodynamic properties can then be derived 
from r2 via the standard thermodynamic relations. 

In a practical calculation, it is more convenient to focus 
on how the virial coefficients are affected by interactions. 
We then may consider the differences AQ„ = Q„ — Qn 



and A&„ 



r,(l) 



where the superscript "1" denotes 



an ideal, non-interacting system having the same fugac- 
ity. As noted in the previous section, our spectrum of the 
eigenstates does not include the non-interacting solutions 
to the boundary value problem. We deal with this issue 
by removing these states from both the interacting and 
non-interacting summations that make up the trace dif- 
ferences AQn. Since they have the same energy with or 



without interactions, this does not affect our results. Ac- 
cordingly, we may rewrite the thermodynamic potential 
in the form, 

n = r!(i) - keTQi [M2z^ + ■■■ + Ab„z'' + •••], (27) 

where fl^^' is the non-interacting thermodynamic poten- 
tial with the same fugacity and 

A62 = Ag2/Qi, (28) 

A63 = AQ3/Q1-AQ2, etc. (29) 

We now describe how to calculate the non-interacting 
thermodynamic potential 57*-^-' and the virial coefficients 
Abn. 



A. Non-interacting thermodynamic potential fi'^' 

Let us consider a two-component non-interacting 
Fermi gas in the thermodynamic limit. In the limit of 
a large number of fermions, the non-interacting thermo- 
dynamic potential fJ'-^ns given semiclassically by, 



f^w = -4 



dpd]i 



13 J {27rr 



= -2 



{hujf 



In 



tlnfl 



- + f-w-p -/i) 



dt. 



m 

(31) 



Subsequently, the number of atoms, N'^^'^ = —957^^^9/1, 
and the entropy, S'^'^^ = —dfl'^^^/dT, may be calculated, 
as wefi as the total energy, E^^^ = Q^^^ + TS^^^ + /j.N'^^I 
We find that. 



and 



ATd) 



^(1) _o(^sr) 



keT 



ze 



1 + ze- 



-dt 



(32) 



E^'> = 2- 



{hujy 



ze 



1 + ze- 



-dt^-2Vt^^\ 



(33) 



B. Second virial coefficient A62 

We now calculate the second virial coefficient. We are 
interested in the limit of a large number of fermions 
{N ^ 1), a situation that will mostly likely happen 
in experiment. As the Fermi energy Ep or the Fermi 
temperature Tp = Ep/ks is given by Ep = N^/'^fiw 
and the temperature T ^ Tp, we shall define a re- 
duced trapping frequency a) — huj/kBT <C 1. The ther- 
modynamic limit is reached in the limit of w — >■ . 
In this limit, the single-particle partition function, de- 
termined by the single-particle spectrum for a 2D har- 
monic oscillator Enm = (2n -I- \m\ + l)fta; is given by 



Oi = 2/(e+'^/2 - e-"/2) ~ 2 {kBTy / {hwy , which can 
also be determined from the first-order expansion of the 
non-interacting thermodynamic potential $7^^^ . The pref- 
actor of two accounts for the two possible spin states of 
a single fermion. 

The second virial coefficient A62 is given by AQ2- It 
is readily seen that the summation over the center-of- 
mass energy in Q2 gives exactly Qi/2. Using the relative 
two-body energy Erei = (2i^„ + l)huj, where i^„ is n-th 
solution of Eq. ([T]), we find that, 



A&2 



1 



Y^ 



2^^ 



-(2i/„ + 1)lD _ -(2i/<,i' + l)<i 



,(1) 



(34) 



where the non-interacting vk — n in = 0,1,2,...) is a 
non-negative integer. 



Third virial coefficient Afos 




Figure 4: (Color online). Second and third virial coefficients 
as a function of the interaction strength EbJEf at different 
temperatures, T/Tf ~ 0.5 (solid lines), 1.0 (dashed lines), 
and 2.0 (dot-dashed lines). 



The third virial coefficient, given by A63 = AQ3/Q1 — 
AQ2, is more difficult to calculate. Both the term 
AQ3/Q1 and AQ2 diverge as w — ^ 0, with the lead- 
ing divergences canceling each other. We thus have to 
separate out carefully the leading terms and treat them 
analytically. It is easy to see that the spin configura- 
tions of tit a-iid Iti contribute equally to Q3. As Qi 
in the denominators cancels exactly with the summa- 
tion over the center-of-mass energy, we have AQ3/Q1 = 
[J2eM~Erei/kBT)-J2eM-El.]}i/kBT)]. To calculate 
this, it turns out to be important to analyze the behavior 
of Erei at large energies. 

To this aim. we define a relative energy Erei, which is 
the solution of Eq. p^ without the exchange term Cnm- 
The utility of Erei is that it can be constructed directly 
from the two-body relative energy. In the subspace with 
a total relative momentum m, it takes the form 



Erei = (2n +\m\ + l)hio + (2i/ + l)hu!, 



(35) 



where v is the solution of the two-particle spectrum Eq. 
([7]) . At large energies where the exchange effect becomes 
less important, the full spectrum Erei approaches Erei 
asymptotically. There is an exception, however, at zero 
total relative momentum m — 0. The solution of Erei at 
n = and to = is spurious, due to the exchange op- 
erator which leads to a vanishing relative wave function. 
It therefore cannot match any solution of Erei- In the 
m = subspace, we must require n > 1 in Eq. (I35p . 

Interestingly, if we retain the spurious solution in the 
m = subspace, the difference [^exp(— E'rei/fesT) — 
^exp(— E'^gj/fc^T)] gives AQ2 exactly, since the first 
part in Eq. (|35p is identical to the spectrum of center-of- 
mass motion. The spurious solution gives the contribu- 
tion. 



n 



-(2t/„+2)w _ -(2iy^i'+2)LD 



which should be subtracted. We thus finally arrive at the 
following expression for the third virial coefficient. 



A6, = 



n 



^-E^^i/ksT _ g-Er^i/kgT 



2e-^Ab2. (37) 



2e-^Ab2, (36) 



The summation should be taken over all the possible rel- 
ative energy levels Erei and their asymptotic counter- 
parts Erei- It is well-behaved at arbitrary interaction 
strengths. 



D. Numerical results of virial coefficients 

We have numerically calculated the second and third 
virial coefficients as functions of interaction strength and 
temperature, with a small reduced trapping frequency 
w <C 1. To ensure the accuracy of the calculations for 
A63, we typically use a hundred thousand relative ener- 
gies Erei - The dependence of the virial coefficients on uj 
may be removed by a careful scaling analysis. 

Fig. 4 shows the evolution of the virial coefficients 
with increasing interaction strength, as characterized by 
the dimensionless two-body binding energy Eb/ Ep. The 
coefficients diverge exponentially in the strongly attrac- 
tively interacting limit, due to the formation of tightly 
bound molecules. The lower the temperature, the faster 
the divergence. 

Fig. 5 presents the temperature dependence of the 
virial coefficients at two interaction strengths, Eb = 
Q-2Ep and Eb = 0-lEp. The coefficients vary strongly 
with the temperature in the degenerate regime (T < Tp). 
However, approaching the high-temperature Boltzmann 
limit (T ^ Tp), the coefficients tend to saturate to a 
semiclassical value. 




Figure 5: (Color online). Temperature dependence of the sec- 
ond and third virial coefficients at two interaction strengths, 
Eb = 0.2Ef (soUd lines) and Eb = O.lEp (dashed lines). 



IV. HIGH-T THERMODYNAMICS OF 
STRONGLY CORRELATED FERMIONS IN 2D 

We are now in position to study the equation of state in 
the high temperature regime. Using the thermodynamic 
relations, it is easy to obtain. 



and 

E= -2n 






[2Ab2z^ + 3Ab3z^ + ■ ■ ■] , (38) 



, {kBTf T 
' {hujf 



Tp 



[Ab', 



z'= + . 



(39) 



where we have defined A6Jj = d(Abn)/d{T/Tp) and the 
Fermi temperature Tp = yNhuj/kB- The entropy is 
then calculated by using S — {E ~ il — ^N)/T, where 
H = ksTlnz. Eqs. ^, ^, and §U^, together with 
the non-interacting number equation (|32|) , form a closed 
set of expressions for thermodynamics. 

We perform the calculation at a given fugacity within 
the trap units h — m = uj — ks — 1- In the case of 
thermodynamic limit, the temperature is fixed to an ar- 
bitrary constant (i.e., T — 100). The virial coefficients 
and their derivative with respect to the reduced temper- 
ature are known as the input. We then calculate N by 
using the number equation pSI) with an initial guess of 
the reduced temperature T /Tp and obtain in turn the 
Fermi temperature Tp = ^/N . The reduced temperature 
T/Tp is updated. We iterate this procedure until the 
final number of fermions and the reduced temperature 
converges within a given relative error. We then calcu- 
late the total energy using Eq. ((5^ and consequently 
the entropy S = {E - fl)/T - Nks Inz. We finally plot 
the chemical potential, entropy or energy per particle, 
pijEp, S/{NkB), and E/{NEp), as a function of the re- 
duced temperature T/Tp. 
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Figure 6: (Color online). Temperature dependence of the 
chemical potential, entropy, and energy of a strongly corre- 
lated Fermi gas in a 2D harmonic trap. The predictions of 
virial expansion up to the third- and second-order are shown, 
respectively, by the solid lines and dashed lines. For com- 
parison, we also show the ideal, non-interacting results using 
dot-dashed lines. 



Fig. 6 gives the high-temperature equations of state 
of a strongly correlated 2D Fermi gas at a typical inter- 
action strength Eb = 0.2Ep. Compared with the ideal, 
non-interacting results, the equations of state of a 2D 
trapped Fermi gas are strongly affected by interactions, 
even in the high temperature regime. The applicability of 
the quantum virial expansion method may be examined 
by comparing the prediction of expansions of different 
orders. We estimate conservatively that the third-order 
virial expansion is reliable down to the Fermi degeneracy 
temperature, T ~ Tp. 



V. CONCLUSIONS AND REMARKS 

In conclusion, we have presented the exact three- 
particle energy eigenstates in a two-dimensional har- 
monic trap, for identical interacting fermions and bosons. 
The energy spectra have been discussed in detail. We 
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have identified two types of energy levels, one contain- 
ing a molecule and the other consisting of individual 
atoms. The latter branch may be interpreted as the en- 
ergy spectrum of a repulsively interacting system. For 
three strongly interacting bosons, we have found two uni- 
versal three-body bound states, corresponding to a self- 
bound boson droplet. The calculated binding energy of 
the droplet is in reasonable agreement with a previous 
theoretical prediction |35| . 

Based on the these exact solutions, we are able to pre- 
dict the high-temperature thermodynamics of a strongly 
correlated quantum gas, by applying a quantum virial 
expansion method. We have calculated for the first time 
the second and third virial coefficients of a strongly corre- 
lated two-dimensional Fermi gas in a harmonic trap and 
have calculated in turn the temperature dependence of 
the chemical potential, entropy and energy. Motivated by 
the striking experimental confirmation of quantum virial 
expansion prediction for strongly interacting fermions in 
three dimensions [50|, we anticipate that our prediction 
in two dimensions will be tested in future experiments 
of two-dimensional Fermi gases. Our thermodynamic 
results may also provide a useful benchmark for future 
quantum Monte Carlo simulations at high temperatures 
for two-dimensional systems of ultra-cold atoms. 



Here, for convenience we have set d ~ 1 as the unit 
of length. Ln is the generalized Laguerre polynomial 
and U is the second Kummer confluent hypergeometric 
function. A direct integration for C„„' is difficult, since 
the second Kummer function becomes singular close to 
the origin. Moreover, the integration for different val- 
ues of Vm.n' makes the numerical calculation very time- 
consuming. 

Thus, it is better to use a different strategy by writing. 







Here, we have used the mathematical identity. 



(A4) 



T(-v)U{~v,l,x^) 



v^ ^k 



L,. ix^) 



^-^ k — v 

k=0 



(A5) 



Therefore, we arrive at 
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L-n 



'^ k- 



fc=0 



V2 






'fc7 



where 



(A6) 



Appendix A: Calculation of Cnn' 

In this Appendix, we outline the details of how to con- 
struct the matrix element C„„' in Eq. pU)) . which is 
given by. 



n 



pdpR„ 



(P) Rn',n (2) 



R 



'V3 



fcO 



(A7) 



V3 



C„ri' = / pdpRnm (p) Rn'm f ^ ) '^2p {-^ Pi '^rn,n') , 



where 



Rnm [p) 



2nl 



[n + |m|) 



can be calculated with high accuracy by using an ap- 
propriate integration algorithm. We note that, with a 
cut-off TT-max for the number of expansion functions (i.e., 
n, n' < rimax), C*^,^ vanishes identically for a sufficient 
large k > fcmax ~ 4nniax- Thus, the summation over k in 
(Al) Eq- (|A6p terminates naturally and one does not need to 
worry about the convergence problem. 

In the practical calculation, we tabulate and store the 
pl"i|g-p Z'^jJiiA fp2\ ^ (A2) coefficients C^/j, in a file, for some given total relative 

angular momentum m. Thus, the calculation of C„„' for 
different values of i'rn,n' reduces to a simple summation, 
which is very efficient and fast. We confirmed numerically 
that the matrix C„„' is symmetric, i.e., C„„' = C„'„. 
A standard diagonalization algorithm can therefore be 



is the radial wave function of an isotropic 2D harmonic 
oscillator and the two-body relative wave function 



V'2p — r(— ;/„„')[/(— i^„j„', 1, -p ) exp( — p~). (A3) adopted for the matrix A/ or Ab. 
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